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CONVECTION WITH NO-SLIP BOUNDARY 
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. CRITICAL LENGTH SCALES 
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^ ' Abstract. We study the Rayleigh-Benard convection in a 2-D 

pL^ . rectangular domain with no-shp boundary conditions for the ve- 

locity. The main mathematical challenge is due to the no-slip 
boundary conditions, since the separation of variables for the linear 
eigenvalue problem which works in the free-slip case is no longer 
possible. It is well known that as the Rayleigh number crosses a 
' critical threshold Rc, the system bifurcates to an attractor, which 

is an (m— 1) -dimensional sphere, where m is the number of eigen- 
values which cross zero as R crosses i?^. The main objective of 
this article is to derive a full classification of the structure of this 
bifurcated attractor when m = 2. More precisely, we rigorously 
prove that when m = 2, the bifurcated attractor is homeomor- 
phic to a one-dimensional circle consisting of exactly four or eight 
1^ , steady states and their connecting heteroclinic orbits. In addition, 

' we show that the mixed modes can be stable steady states for small 

, Prandtl numbers. 
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1. Introduction 

The Rayleigh-Benard convection problem is one of the fundamental 
problems in the physics of fluids. The basic phenomena of the Rayleigh- 
■ Benard convection in horizontally extended systems are widely known. 

The influence of the side walls, although not studied as throughly as the 
horizontally extended case, is of practical importance for engineering 
applications. 

In this paper we study the Rayleigh-Benard convection in a 2-D 
rectangular domain with no-slip boundary conditions for the velocity. 
This problem is also closely related to the problem of infinite channel 
with rectangular cross-section which has been studied by Davies- Jones 
|DJ70j . Luijkx-Platten |LP81] and Kato-Fujimura |KF00j among oth- 
ers. 

The linear aspects of the problem we consider in this paper have 
been studied by Lee-Schultz-Boyd |LSB89j . Mizushima |Miz95] and 
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Gelfgat |Gel99] . In these papers, the critical Rayleigh number and the 
structure of the critical eigenmodes have been studied for small aspect 
ratio containers. 

From dynamical transition and pattern formation point of view, Ma 
and Wang |MW04l IMW07] proved that under some general bound- 
ary conditions, the problem always undergoes a dynamic transition to 
an attractor S^j as the Rayleigh number R crosses the first critical 
Rayleigh number Rc. They also proved that the bifurcated attractor, 
homological to S"^~^, where m is the number of critical eigenmodes. 

In the 2-D setting that we consider, m is either 1 or 2 and the latter 
case can only happen at the critical length scales where two modes with 
wave numbers k and k+1 become critical simultaneously. When m = 1, 
the structure of E/j is trivial which is merely a disjoint union of two 
attracting steady states. Thus our task in this paper is to classify the 
structure of the attractor when m = 2. This has been studied recently 
in |SW13j for the 3D Rayleigh-Benard problem where the boundaries 
were assumed to be free-slip for the velocity and the wave numbers of 
the critical modes were assumed to be equal. 

The main mathematical challenge in this paper is due to the no- 
slip boundary conditions since the separation of variables for the linear 
eigenvalue problem which works in the free-slip case is not possible any- 
more. To overcome this difficulty, the main approach for our study is 
to combine rigorous analysis and numerical computation using spectral 
method. 

As we know, spectral methods have long been used to address the hy- 
drodynamic instability problems. In fact, in his seminal work |0rs71j . 
Orszag studied the classical Orr-Sommerfeld linear instability problem 
using a Chebyshev-tau method. In this paper, to treat the linear eigen- 
value problem, we employ a Legendre-Galerkin method where compact 
combinations of Legendre polynomials, called generalized Jacobi poly- 
nomials, satisfying all the boundary conditions are used as trial func- 
tions. The main advantage of our Legendre-Galerkin method is that 
the resulting matrices are sparse which allows a very efficient and accu- 
rate solution of the linearized problem; see also Hill-Straughan |HS06j 
and Gheorghiu-Dragomirescu |GD09j . 

Once the eigenpairs of the linear problem are identified, the transi- 
tion analysis is carried out by reducing the infinite dimensional system 
to the center manifold in the first two critical eigendirections. The co- 
efficients of this reduced system are calculated numerically. Our main 
results are described below. 

We first classify the eigenmodes into four classes according to their 
parities using the symmetry of the problem. Then we numerically show 
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that the first two unstable modes are always parity class one or two. 
Then we study the transition near the critical length scales where two 
eigenvalues become positive simultaneously. Next, we rigorously prove 
that the local attractor at small supercritical Rayleigh numbers is in 
fact homeomorphic to the circle which has four or eight steady states 
with half of them as stable points and the rest as saddle points. The 
critical eigenmodes are always bifurcated steady states on the attractor 
and when the attractor has eight steady states, the mixed modes which 
are superpositions of the critical eigenmodes are also bifurcated. 

Second, let /3i and /32 denote the two largest eigenvalues of the lin- 
earized problem. We find that a small neighborhood of /3i = = 
in the /3i-/32 plane can be separated into several sectors with different 
asymptotical structures. In particular, we find that there is a critical 
Prandtl number Pre for the first two critical length scales L = 1.5702 
and L = 2.6611, such that for Pr < Pre, there is a sector in this plane 
for which mixed modes are stable fixed points of the attractor. For 
Pr > Pre, the mixed modes are never stable and instead there is a 
sector in this plane in which both of the critical eigenmodes coexist as 
stable steady states. In this case, the initial conditions determine which 
one of these eigenmodes will be realized. The critical Prandtl number 
is around 0.14 for the first critical length scale L = 1.5702 and around 
0.05 for the second critical length scale L = 2.6611. For higher critical 
length scales we found that mixed modes are never stable points of the 
attractor. 

Third, recently Ma- Wang has developed the dynamic transition the- 
ory to study transition and bifurcation problems in nonlinear sciences; 
see [MWj . This paper is a first attempt to combine this theory with the 
numerical tools of the spectral methods to study the detailed structure 
of the transition and pattern formation. 

The paper is organized as follows. In Section 2, the governing equa- 
tions and the functional setting of the problem is discussed. In Section 
3, linear eigenvalue problem is studied. Section 4 states the main theo- 
rem. Section 5 is devoted to the proof of the main theorem. In Section 
6, we demonstrate a method to compute the coefficients of the reduced 
system. And the last section discusses the results obtained by our 
analysis. 

2. Governing Equations and the Functional Setting 

Two dimensional thermal convection with no-slip, perfectly conduct- 
ing boundaries can be modeled by the Boussinesq equations. The gov- 
erning equations on the rectangular domain Q = (0, L) x (0, 1) G 
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read as 

^ + (u ■ V)u = -Pr(Vp + Au) + VRVP^Ok, 
V-u = 0. 

Here u = [u, w) is the velocity field, 6 is the temperature field and p 
is the pressure field. These fields represent a perturbation around the 
motionless state with a linear temperature profile. The dimensionless 
numbers are the Prandtl number Pr and the Rayleigh number R which 
is also the control parameter, k represents the unit vector in the z- 
direction. 

The equations ([T]) are supplemented with no-slip boundary condi- 
tions for the velocity and perfectly conducting boundary conditions for 
the temperature. 

(2) u = ^ = 0, on dn. 

For the functional setting, we define the relevant function spaces: 
H = {{u,9) G (fi,M3) . v - u = 0,u-r2 \9n= O} , 
Hi = { (u, 6) G {Q, M=^) : V ■ u = 0, u |sn= 0, 6 {90= O} . 

For = (u, 9), let G : Hi ^ H and Lpt : Hi ^ H he defined by 

= (^P(PrAu + /ffVP^^k), \^Vf^w + A^) , 

G(0) = -(P(u-V)u, (u-V)^)), 

with V denoting the Leray projection onto the divergence-free vectors. 

The equations ([I])-© supplemented with initial conditions can be 
put into the following abstract ordinary differential equation: 

(5) ^ = L^0 + G(0), 0(0) = 00. 

For results concerning the existence and uniqueness of solutions of 
([5]), we refer to Foias, Manley, and Temam |FMT87] . 

Finally for 0j = (ui,^j), Uj = {ui,Wi), i = 1,2,3 we define the 
following trilinear forms. 

(^(01, 02, 03) = - / (ui ■ V)u2 ■ u^dxdz - / (ui ■ V)6'2 ■ Osdxdz, 

(6) Jn Jn 

Gs{4>l, 02, 03) = ^(01, 02, 03) + ^(02, 01, 03)- 
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3. Linear Analysis 
We first study the eigenvalue problem Lr</) = /30 wliicli reads as 

Pr(Aw - ^) = 

Pr(Au7 - ^) + VrVf^O = (3w, 

(7) 

divu = 0, 

u = e = at dVt. 

Below we list some of the properties of this eigenvalue problem. 

1) The linear operator is symmetric. Hence the eigenvalues /3j 
are real and the eigenfunctions 0j are orthogonal with respect 
to L^-inner product. Moreover there is a sequence 

< i?i < i?2 < ■ ■ ■ 

such that (3i{Ri) = 0. Ri is found by setting /3 = in ([7]). In 
this case the problem becomes an eigenvalue problem with a/R 
as the eigenvalue. 

2) We have 

(3i{R) =OiiR = Ri, 

which can be seen by computing the derivative of Pi with respect 
to i? at i? = i?,-. 



(9) 



where {ui,Wi,9i) is the ith eigenfunction. Also at i? = Ri, by 
the third equation in (0), u;^ = -R~^^^Fi~^/^Aei as l3i{Ri) = 0. 
Plugging these into (IHl) and integrating by parts, we see that 



dPi 

I R=Ri 



In 


iV^il 


"^dx 


R^jui\u^\ 


|2 + 


\Wi 


|2+ |^i|2)c/x 



> 



dR 

3) We denote the critical Rayleigh number R^ = Ri. That is 

{< ifR<Rc, 

= iiR = Rc, i = l,...,m 

> iiR> Rc. 
/3i{Rc) < 0, i>m. 
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Class 1 


Class 2 


Class 3 


Class 4 


^(e,e), 9{o,e) 


^(o,e), 9{e,e) 


^(e,o), 9{o,o) 


^^(0,0), 61(6,0) 



Table 1. Possible parity classes of the eigenf unctions 
of the linear operator. 



m in ([9]) does not depend on the Prandtl number Pr but only on 
L. To see this, one makes the change of variable 9 = \/Pt9' so 
that the solution of ([7]) for the eigenvalue of /3 = is indepen- 
dent of Pr. By simplicity of the first eigenvalue (see Theorem 
3.7 in Ma- Wang |MW05] ). for almost every value of L except 
a discrete set of values, m in ([9]) is 1. 

Introducing the streamfunction tpz = u, tpx = —w, we can eliminate 
the pressure p from the linear eigenvalue problem ([7]). 



PrA^V; - VrVf^9^ = (3{R) Aip, 
(10) -y/nVP^^, + A9 = f3{R)9, 

ijj = ^ = 9 = on dn. 
on 

The linear equations ([7]) satisfy several discrete symmetries which 
may be found from the known non-trivial groups of continuous Lie 
symmetries of the field equations ([1]); see (Hydon |HydOO| , Marques- 
Lopez-Blackburn |MLB04] ). However, for the problem we consider, it 
can be easily verified that the linear equations have reflection symme- 
tries about the horizontal and vertical mid-planes of the domain. Thus 
we can classify the solutions of the linear problem into four classes 
with different parities which are as defined in Table [T] where, for ex- 
ample, ip{o, e) means that ip is odd in the x-direction and even in the 
^-direction. 

We will employ a Legendre-Galerkin method (cf. Shen |She94] . 
Shen-Tang-Wang jSTWllj ) to solve the linear eigenvalue problem 
f lTOj) . For this, first we transform the domain with the change of vari- 
able 

2x 

(x, z) G (0, L) X (0, 1) ^ (X, Z) = i—- 1, 2z - 1) G (-1, 1)1 

The approximate solutions {tp^ , 9^) of ffTOl) will be sought in the finite 
dimensional space 

= span{{ej{x)ek{z), fi{x)fm{z)) \ j, k,l,m = 0,1, . . . , N - 1}, 
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where Cj and fj are generalized Jacobi polynomials (cf. Guo-Shen- 
Wang |GSW06] . Shen-Tang-Wang |STW11] ) which satisfy the bound- 
ary conditions 

e,(±l) = De,(±l) = /,(±l)=0. 

Here D denotes the derivative. The polynomials and are defined 
as in Chapter 6 of Shen-Tang-Wang |STWllj . 

(11) h{z) = U{z)-U+^{z), 

(12) e\z) = — ^ ~ ^»+2(^) + ^^i+ii^) 



(2z + 3) (4z + 10)1/2 

where Lj is the i^^ Legendre polynomial. The coefficient of Cj guarantees 
that {D'^ei,D'^ej) = 6ij. 

We write the approximate solutions of the equation (fTO!) with coef- 
ficients to be determined by 

(13) = E E ^fke,{x)e,{z), = ^ E ^U(^)fki^)- 

Here N = 2NxNz denotes the total degrees of freedom. 
Let us define for i, j = . . . , m — 1, 

(A^),, = (I)'e„ D^e,) = {A',^),, = {D^e„ e,) = -{De„ De,), 

(^3 )«i ~ i^iy ^i)' (^4 )ij — i^iy fj)y 

and for j = 0, . . . , A^^. - 1, A; = 0, . . . , A^^ - 1, 

Using the following property of the Legendre polynomials 

(14) (2z + 3)L,+i = D(L,+2-i^.), 

it is easy to see that: 
(15) 



Li+3 — Li^i 2 / \ / 2i + 5 



De,{z) = D'e,{z) = J ^—L,^,{z), Df, = -(2z+3)L,+i. 

V4« + 10 V i 

By (fT5|) . it is easy to determine the elements of the matrices Ak by the 
orthogonality of the Legendre polynomials. In particular, the matrices 
A^, . . . A™ are banded, and except for A^ and A!^, they are symmetric. 
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Putting (fT3|) into ( fTOl) . multiplying the resulting equations by em{x)en{z) 
and fm{x)fn{z) respectively and integrating over — l<x<l,— 1< 
2; < 1 to obtain 



(16) 

Here: 
(17) 



B 



N 



N 



D 



where 



PrXi 
X3 

X4 
X, 



N 



NxN 



c 



PrX^^ 



PrX2 




NxN 



NxN 



N\ 



1^ 

^4 3 



vec 



J iVxl 



aN^ 

^3 ) 



X2 = Y^A ® (^7 

22 



X, = -,A, 



N, 



v4f^ + 2' A 



2 aN, 

5 



A 



N^ 
6 5 



X4 = ^A^^ ® 



2 +2M^ 



^3 ' 



X. = A 



N, 



aN^ 
^6 



In (1171) and f|T8|) we use the following notations. For a m x k matrix 
M, vec(M) is the mk x 1 column vector obtained by concatenating the 
columns M, of M, i.e. 



vec 



Ml M2 



) 



Ml M2 



r 



stands for the zero matrix and A® B = {ajji?}jj=o,i, ...,(? is the Kro- 
necker product of A and B. To obtain (flGll . we used the following 
properties of the Kronecker product. 

vec{AXB) = (fi^ ® y4)vec(X), {A ® 5)^ = ® 5^. 

We note that the matrices B^ , and in (fT6|) are sparse, 5^ 
and are symmetric while is skew-symmetric. 
From our linear analysis, we find the following results. 

• Our numerical analysis suggest that = 6 + 2k ^ 6 + 2L and 
A^z = 8 is enough to resolve the critical Rayleigh number and the 
first critical mode which has k rolls in its stream function. We 
have checked that increasing and Nz by two only modifies 
the fourth or fifth significant digit of the result. 
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Figure 1. 

if) (on top) and 9 (on bottom) of the first critical mode on 
(-L,L) X (-1,1) for L = 1,...,5. 




Figure 2. 

The left figure shows the plot of the first critical stream function for 
L = 1 (the top left plot in Figured]). The right figure shows the 
enlarged plot at the corner. 

• In Figure [H the first critical mode is shown for the length scales 
L = 1, . . . , 5. Note that the first critical stream function and 
the temperature distribution has always even parity in the z- 
direction while their x-parity alternates between odd and even 
as the length scale increases. As observed in Mizushima |Miz95] , 
we also verify the existence of the Moffatt vortices on the corners 
of the domain which are due to corner singularities as shown in 
Figure [21 

• For L < 21 we observed that m in ([2]) is either 1 or 2. More- 
over, m = 2 only at the critical length scales which are given 
in Table [2l The results found are in agreement with those in 
Mizushima |Miz95j and Lee-Schuhz-Boyd |LSB89j . 

• The marginal stability curves of the first few critical eigenvalues 
are given in Figure [31 The figure demonstrates that the parities 
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T 

w 


k 


Kc 


"I\.T 


T 


k 


1 ) 

tXc 




1.5702 


1 


3086.6554 


8 


6.7711 


6 


1764.3754 


18 


2.6611 


2 


2113.776 


10 


8.7992 


8 


1740.9174 


22 


3.7048 


3 


1906.3395 


12 


10.8229 


10 


1729.5398 


26 


4.7329 


4 


1826.4099 


14 


15.8738 


15 


1717.805 


36 


5.7541 


5 


1786.8833 


16 


20.9197 


20 


1713.5226 


46 



Table 2. At L = Lc, two modes become unstable. One 
of the modes has k and the other one has k + 1 rolls in x- 
direction in their stream functions. The critical Rayleigh 
number at this length scale is Re- and = '8: are 
the number of polynomials used in the x and z directions 
respectively. 




Figure 3. The marginal stability curves of the first few 
eigenvalues with eigenf unctions of different parity classes. 



of the first two critical modes can only be of parity class one or 
two as given by Table [H 

Also it is seen in these figures that there is a repulsion of 
the eigenvalues. Namely the neutral stability curves of the 
same parity type do not intersect each other. Such a repulsion 
does not occur for free-slip boundary conditions. This repulsion 
arises from a structural instability of the transform of matrices 
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into a Jordan canonical form and a detailed analysis can be 
found in Mizushima-Nakamura |MN02j . 

4. Main Theorem 

Let m be the number of modes which become critical as the first 
Rayleigh number Rc is crossed as given by ([9]). Ma and Wang |MW04t 
IMW07] proved that under some general boundary conditions, the prob- 
lem has an attractor S/j which bifurcates from (0, -Rc) as R crosses 
Rc. They also proved that the dimension of the bifurcated attractor is 
m — 1 < dim(S/j) < m. When m = 1, the structure of is trivial 
which is merely a disjoint union of two attracting steady states. 

As stated before, in our problem m is either 1 or 2. And the latter 
case can only happen at a critical length scale Lc where two eigenmodes 
with consecutive wave numbers become critical. 

Numerically, it turns out that the critical Rayleigh numbers for 
modes with parity 3 or 4 are much greater than those for modes with 
parity 1 or 2. This can be seen from the Figure [31 

We will assume the following. 

' 1. (/3i,0i), (/32,02) are the first two critical eigenpairs. 
2. 01 has wave number k, 02 has wave number k + 1 
(19) < where k is a positive integer. 

]. One of the eigenmodes {0i, 02} is of parity class 1, 
while the other is of parity class 2 as given in Table [T] 

Let Hi and y2 be the amplitudes of 0i and 02 respectively. Then 
in the proof of the main theorem, we show that the dynamics of the 
system close to R = Rc and L = Lc is governed by the equations 

= f^iVi + yiiauvl + flisl/a) + o(3), 

(20) 

= /321/2 + 1/2(022?/? + a24yl) + 0(3). 

Here /3j is the eigenvalue corresponding to mode i, and o(3) denotes 
o(3) = o(|(yi,2/2)|') + 0(|(2/i,y2)|'max|/3(i?)|). 

4 = 1,2 

Let us define 

(21) -Di = a22/3i-aii/32, ^2 = 013/^2-024/31, -03 = 011024-013022- 

To state our main theorems, we assume the following non-degeneracy 
conditions 

(22) an ^ 0, 024 ^ 0, ^ 0, D2 ^0,0^^ 0. 
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(a) Di < 0, L>2 < (b) L»i > 0, L»2 > (c) Z^i < 0, L>2 > (d) Di > 0, L»2 < 
0, Da < 0, D3>0 



Figure 4. Transition Scenarios. 

Finally, let us define the following. 
(23) 

ifi = (— i = 1,2, (modes with wavenumber k) 

tpi = (— 1)*02, = 3,4, (modes with wavenumber k + 1) 

ipi = Cj0i + (ii</'2, i = 5,6, 7, 8, (mixed modes) 

where = Cq = —Cj = — Cg and = —d^ = d^ = —dg. 

Theorem 4.1. Under the assumptions (fT9|) and ( 122|) . there is an at- 
tractor bifurcating as R crosses Rc which is homeomorphic to the 
circle when L is sufficiently close to a critical length scale L^. More- 
over, consists of steady states and their connecting heteroclinic or- 
bits. Let n{T,fi) denote the number of steady states on S/j, S denote the 
stable steady states andU denote the unstable steady states on upto 
topological equivalency. Then we have the following characterization of 
which is also illustrated in Figure^ 



(i) IfD,<0, D2<0, Ds<0, 


^^(Sfl) = 


8, S 


= {(^i\i = l. 


2,3,4} 


U = {if, 1 z = 5,6,7,8}. 










(ii) IfD,>0, D2>0, 1^3 >0, 


'^(Sr) = 


8, S 


= {ipi\i = 5, 


6,7,8} 


U = {ip,\i = 1,2,3,4}. 










iii) IfD, <0, D2> 0, ^(Sr) 


= A,S = 


{Vi 


1 i = 1,2}, U 


= {Vt 


z = 3,4}. 










iv) IfDi >0, D2< 0, ^(Sk) 


= A,S = 




1 i = 3,4}, U 


= {Vt 


z = l,2}. 











According to Theorem 14. 1^ the structure of the attractor depends 
on Di, D2 and -D3 which in turn depends on the coefficients of the 
reduced equations. By (12T]) . D3 has a definite sign whereas Di and D2 
vanish at the criticality Pi = (32 = 0. In the proof of Theorem 14. H we 
analytically prove that the coefficients an, and 024 are negative. Our 
numerical computations indicate that ais is also always negative. We 
observed that 022 and D3 can be both positive and negative. 
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Mixed 
modes 
stable 



First 

eigenmode 
stable 



(a) 022 > 



Second 
eigenmode 



Second ^2>< ^ 
eigenmode 
stable 




Second 

eigenmode 

stable 



D=0 

Second 1 V 
eigenmode / ^"'"'^ 
stable / eig<:nmodes 
" stable 




(b) a22 < 0, Da > (c) 022 < 0, < 



Figure 5. The transition scenarios in /3i-/32 plane. In 
the above cases we assume that ais < which is due 
to our numerical observations. The arrows on the lines 
Di = 0, D2 = indicate in which directions Di and D2 
increase. First and second eigenmodes correspond to the 
eigenmodes with wavenumber k and k + 1 respectively. 

That gives three possible cases depending on the signs of 022 and 
D3. In Figure |5l we classify these cases in a small neighborhood of 
/3i = /32 = in the /9i-/32 plane according to our main theorem and the 
following observations. 

• If 022 > then > 0, Z^i > but D2 changes sign in the first 
quadrant. 

• If 022 < and -D3 > then Di and D2 changes sign in the first 
quadrant. Moreover the case where both Di < and -D2 < is 
not possible. 

• If 022 < and D3 < then again Di and D2 changes sign in 
the first quadrant. This time the case where both Di > and 
-D2 > is not possible. 



5. Proof of the Main Theorem 

We will give the proof in several steps. 
STEP 1. The reduced equations. When there are two critical 
modes 0i, 02? the center manifold is a two dimensional manifold embed- 
ded in the infinite dimensional space. We denote the center manifold 
function by: 

$ = yl^i + yiy2^2 + yl^3 + o{y^), = (^'i, 9^)^ . 

To study the dynamics on the center manifold, we plug in 



(24) = + 1/202 + 2/1^1+ yil/2$2 + 1/2 $3 + 0(2) , 
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into ([5]), take the inner product with 0i, 02 and use the orthogonahty of 
the eigenvectors, thanks to the self-adjointness of the hnear operator. 
The reduced equations read 

% „ . 1 
dt 



(25) 



-(G(0),0,), 1 = 1,2 



We normahze the first two eigenfunctions so that 

(01, 0l) = (02,02) = 1- 

Now if we expand the nonhnear terms in ( l25l) . we get 
dyi 

(26) 

dt 
where 

flfcl = ^^(01, $1, 0fc), ak2 = Gs(01, $2, 0fc) + Gs(02, $1, 0fc) 



/3iZ/i + (aii?/i + ai2?/ii/2 + ai3?/i?/2 + ouZ/s) + o(3), 

/32y2 + (021^1 + 022^12/2 + a23yiy2 + «24y2) + ^(3), 



(27) 



afc4 = ^^(02, ^3, 0fc), flfcS = ^^(01, $3, 0fc) + Gs 



$9 



STEP 2. Parities of the center manifold functions. To com- 
pute the center manifold approximation, we will use the the following 
formula which was introduced in Ma-Wang |MW05] . 

-£jj$i = P2G(0i,0i), 

-£r$2 = P2[G'(0i,02)+G(02,0i)], 
-£jj$3 = P2G(02,02). 



(28) 
Here 



LetX 



P2 '■ H ^ E2, Cr = Lr\e2 '■ E2 -> £"2, 

El = span{0i, 02}, E2 = E^. 
{/ G C{n) I /(-a;, z) = ±f{x, z) and /(x, -z) = ±/(x, z)} 



and let s : X — 7- {±1}^ denote the parity function: 

s(/) = (^x(/),5.(/)), 

where 

s,(/) = ±l if/(-x,^) = ±/(x,^), s,(/) = ±l if /(x,-^) = ±/(x,^) 
Let us define for 0j = {ui,Wi,9i), i = 1, 2 the following. 



(29) G'(0„0,) 





(pi,' 


Pj) 




duj 


dui 


^2( 


(pi,' 


Pj) 




dwj 


dwi 


93{ 


(pi,' 


P3) 




^« dx 


dOi 
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The following lemma can be proved using the basic properties of 
parities. 

Lemma 5.1. If (pi = {ui,Vi,9i) G H Hi, i = 1,2 then for i, j, k = 
1,2, 

1) -sigi{(f)h(l)j)) = s{g2{(f)i, = s{g3{(f)i, = {s^{wiWj),-s^{wi% 

2) s{gk{(l)i,(j)j)) = s{gk{(l)j,(l)i)). 

Hereafter without loss of generality we will assume 

(30) 01 is of parity class 1 and 02 is of parity class 2, 

which are as given in Table [1] 

Using the Lemma I5.H we can prove 

Lemma 5.2. Under the assumption (!30l) . 

■5(5'2(01, 0l)) = 5(5-2(02, 02)) = (1, "1), 5(52(01, 02)) = ("1, "l)- 

Lemma 5.3. Under the assumption fl30|) . P2G{(f)i,(j)j) = G{(j)i,(j)j) for 
^,J = 1,2. 

Proof Note that P2G{(j)i,(j)j) = G{(j)i,(f)j) if (G(0i, 0^), 0^) = for 

i,j, k = 1,2. Now 

(31) 

(G(0i, 0j), 0fc) = / (51 (0i, 4>j)uk + 52(0i, 0i)w^fc + 53(0i, 4>jWk) dxdz 



By Lemma [5.11 and Lemma [5.21 5i(0j,0j) is even in the ^-direction 
while g2{4>i,4>j) and g3{(f)i, are odd in the ^-direction. Since Uk is 
odd and Wk and 6k are even in the 2;-direction, the integral in ([31]) must 
vanish over Q. □ 

Thus by the Lemma [531 and the equation f[28|) . $j = {Ui,Wi,Qi), 
{i = 1,2,3) are solutions of 

-£^$i = G(0i,0i), 

(32) - £/?$2 = G(0i, 02) + ^(02, 0i), 

-£^$3 = G(02,02). 

Using the streamfunction = U, = —W, one can eliminate the 
pressure from these equations to obtain 

PrA ^ - Vi?vPr— = hi := — — + — — , 

ox oz ox 

(33) -^^^ + Ae = h2 := -53, 

ox 

^^J = — = e = on dn. 

on 
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01 


h 


$1 


$2 


$3 


V^i 








^1 


Oi 






^3 


63 


(e,e) 


(o,e) 


(o,e) 


(e,e) 


(o,o) 


(e,o) 


(e,o) 


(o,o) 


(o,o) 


(e,o) 



Table 3. Parities of the first two critical modes and 
the center manifold functions. 



Lemma 5.4. Under the assumption (1301) . the center manifold func- 
tions have the parity as given in TablelBi 



Proof. We can eliminate from the first equation of ( 1331) to obtain 

dh 

PrA^* - i?Pr— ^ = Ahi + ^/rVf^- ^ 



Ae = /i2 + v^Vp^— . 

ox 

Now using Lemma [Sm and Lemma lS^ we see that s{^) = {—Sx{g2), Sz{g2)) 
and s(e) = s{g2). □ 

Using Table [3] we find that the integrands in ai2, ai4, a2i; 0'2z are all 
odd functions of z and hence we have the following. 

Lemma 5.5. Under the assumption (1301) . in (127|) we have 

O12 = Ol4 = '^21 = '^23 = 0. 

As a result of Lemma [5.51 we obtain the reduced equations (l20l) . 

STEP 3. The attractor bifurcation. Now, we will prove that 
the bifurcated attractor is homeomorphic to S^. For this we will need 
the following result. 

Theorem 5.6 (Ma- Wang |MW05] ). Let v he a two dimensional (r > 
1) vector field given by 

(35) vx{x) = f3{X) X — h{x, \), 

forx G M^. Here /3(A) is a continuous function of X satisfying /3{X) = 

for X = Xq and 

< 

h{x,X) = hk{x,X) + o{\x\''), < {hkix,X),x), 

for some odd integer k >3 where hk{-,X) is a k-multilinear field, and 
Ci> Q is some constant. Then the system 

dx/dt = V\{x), X G M, 
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bifurcates from [x, A) = (0, Aq) to an attractor S^ , which is homeomor- 
phic to , for Aq < A < Ao + e, for some e > 0. Moreover, either (i) 
is a periodic orbit, or (ii) consists of an infinite number of singular 
points, or, (Hi) contains at most2{k + l) singular points, consisting 
of 2N saddle points, 2N stable node points and n{< 2{k + 1) — 4A^) 
singular points with index zero. 

Now let 

Kvi^Vi) = {2/1(011^1 + 013^2), 2/2(022^/1+0242/2)}^- 

Lemma 5.7. Assume that $j 7^ for i = 1,2,3. Then for any y = 
(2/1,2/2), 

(36) ih{y),y) = a^y^ + (013 + a22)2/?2/2 + «242/2 < C'bl^, 

where C < 0. 

Proof. 

an = Gs{(j)u 0i) = G(0i, $1, 0i) + G($i, 0i, 0i) 
^ = G(0i, $1, 01) = -G'(0i, 01, $1) = (£^$1, $1). 

Here we used (!32l) and the following properties of Navier-Stokes non- 
linearity 

(38) (0 G(0, 0, 0*) = G(0, 0*, 0), (u) G(0, 0, 0) = 0, 

and —Cr^i = G(0i,0i) which is due fl5^ . 
If we write 

00 

'^'i = X]cj-fc0fc, j = 1,2,3, 

fc=3 

then for j = 1, 2, 3, 

00 

(39) (£^<I>„%) = ^4,/3,||0,,||2<O. 

k=3 

since < for /c > 3 and by assumption there exists A; > 3 such that 
Ci^k 7^ 0. In particular, an < 0. As in flTTI) . we can show that 

^24 = ^,(02, $3, 02) = {Cr<^3, $3) < 0. 

Now if ai3 + a22 < then it is easy to prove fl36|) . Assume otherwise. 
Using fl38|) and fl32l) . we can write 

ai3 = G,(0i, $3, 0i) + G,(02, $2, 0l) 

= ^(01, $3, 0l) + G(<I'3, 01, 0l) + ^,(02, $2, 0l) 
^ ' =-(G(0i),$3)+a(02,<f2,0l) 



(£«$!, $3) +a(02,$ 



18 TAYLAN SENGUL, JIE SHEN, SHOUHONG WANG 

A similar computation shows 

(41) 022 = (/:/j<l>3,$l) +a(01,$2,02). 

Let us define 

(42) a = Gs{(pi, $2, H + Gsih, $2, </>i). 
By (I3H]) and ([32]), 

(43) a = -(G(</>i, <P2) + G((/)2, ^i), $2) = {Cr^2, $2) 

Note that a < by ( !39l) . Using Cauchy-Schwarz inequality and the 
orthogonality of the eigenf unctions, 

00 

(£ij<l>l,$3) = ^/3fcCi,fcC3,fc||0fc||^ 

(44) 



fc=3 

\ 1/2 / 00 \ 1/2 



\k=2, J \fc=3 

= a/ 011024- 

Since, (£r$3, $1) = (£,?$i, $3), we have by gOD-dll, 



ai3 + 0-22 < 2a/Oii024 + 

where a < is given by ( l42l) . Thus, there exists < ei < — on, 
< 62 < — 024 such that 



ai3 + ^22 < 2a/oii024 + « < 2a/ (oii + ei)(o24 + £2)- 
Since 2o6 < o^ + 6^, we have. 



2\/ (an + ei)(o24 + ^2)ylyl < -(^^n + ^1)^^ - (024 + e2)l/2- 
Now, let C = max{— ei, —62}. Then C < and we have 

ih{y),y) < anyf + (013 + a22)yfyj + 0241/2 < ^ix"^ + y'^f. 
That finishes the proof. □ 

Thus by Theorem 15.61 and Lemma [5.71 E^j is homeomorphic to S^. 
Now we will describe the details of its structure by determining the 
bifurcated steady states and their stabilities. 

STEP 4. The steady states and their stabilities. The possible 
equilibrium solutions of the truncated equations of (12(1 follows. 



(45) = (a/^,0), i?2 = (0, J^), M = J^, J^), 

-On V -^24 V -D3 V -D3 



where Di, D2 and -D3 are given by (I2T 
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Due to the invariance of the equations fl2Ul) with respect to {x, y) — )■ 
(— X, y) and {x,y) — > {x,—y), we only consider the positive solutions 
when writing ( 145|) . 

The eigenvalues Ai, A2 of the truncated vector field at the steady 
states Ri, R2 are 

Afi = -2/3i, Af^ = -Di/an, Af^ = -2/32, A^^ = -^2/024. 

Note that Ri is always bifurcated for > 0, i = 1,2. Moreover Ri 
is a stable steady state for /3j > if < for i = 1,2. The trace 
Tr and the determinant Det of the Jacobian matrix of the truncated 
vector field at the mixed states M are 

(46) Tr = ^{a^^Di + 011^2), Det = -^D^D^. 

Notice that the steady states M are bifurcated only when Di, D2, -D3 
have the same sign. Since an and 024 are both negative as shown in 
Lemma 5.7, according to trace-determinant plane analysis, they are 
saddles if Di < 0, D2 < 0, < and are stable if Di > 0, D2 > 0, 

Finally, only the four cases stated in our main theorem can occur. 
To see this, note that according to the Theorem 15.61 and fH5|) - fj46l) . the 
case Di < 0, -D2 < 0, -D3 > is not possible since that would lead to 
only 4 steady states on the attractor which are all stable. Similarly the 
case Di > 0, D2 > 0, D3 < is not possible either which would lead 
to 4 steady states which are all unstable. 

6. Numerical approximation of the coefficients of the 

reduced equations 

To compute the coefficients of the reduced equations fl2U]) . we fix L, 
Pr and R to compute all the eigenvalues and the corresponding 
eigenvectors of ffT^ . 

Numerical computation of the center manifold functions. 
Now we will numerically approximate $1, $2 and $3 which are the 
solutions of the equations (128|) . We will illustrate the method to ap- 
proximate $1 since $2, $3 can be approximated similarly. To determine 
$1, we have to find its stream function \i/ and its temperature function 
G which are determined by equations fl33p . 

Since we do not have hi and /12 in fl33p exactly, we approximate them 
by , h2 as below 

uN J.N qN „i,N qN 

(47) = + = - ^^-^^ 



2) 
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Here {ipi , 6i) is the first critical eigenf unction of the discrete problem 

(48) =J2I1 {^i,mnem{x)er.{z), 9^,^Jm{x)U{z)} . 

m=0 n=0 

The Legedre-Galerkin approximation of the problem ( l33l) . 

As in the linear eigenvalue problem, we discretize the equations ( l33l) 
using the generalized Jacobi polynomials (fTTl) - (fT2l) . 

(49) {vi/^,e^} =J2I1 {^LU^)en{z), eLUx)u{z)] . 

m=0 n=0 

We plug in 6^, /if , /i^ for ^, 6, /ii, /i2 in ([33]) and multiply 
the resulting equations by Jacobi polynomials ej{x)ek{z), fj{x)fk{z) 
and integrate over — l<a;<l,— l<2;<lto reduce fl5^ to the 
following finite dimensional linear equation 

(50) {B^ - \^C^)x = b. 
Here and are given by f[T7|) and 



(51) X 



T 



vec(^^) vec(e^) J^^^^, b= vec(fii) vec(52) 



For < J < A^^. - 1, < A; < A^^ - 1 



(52) 



{Bi)jk = I J ^1 {x, z)ej{x)ek{z)dxdz, 
{B2)jk= I J h2 ix,z)fj{x)fkiz)dxdz. 



Now Cj is a polynomial of degree j ' + 4 and by fHTI) and fHS]) . /if' is 
a polynomial of degree at most (2A^j. + 6,2Nz + 6). Thus the above 
integrands are of degree at most (3A^. + 9, 3A^2 + 9). Since the Legendre- 
Gauss-Lobatto quadrature with A^ + 1 quadrature points is exact for 
polynomials of degree less or equal than 2A^ — 1, the integrals in (l52l) 
can be replaced by the following discrete inner products. 

{Bi)jk= ^ ^ h'^{x^,Zn)ej{x^)ek{zn)u:l,u:n^ 

/rn\ m=0 n=0 

{B2)jk= h2{Xm,Zn)fj{Xm)fk{Zn)uJ^UJn- 
m=0 n=0 
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Here {xj,Wj}j^Q~^^ and {zj,Wj}j^Q^^ are the Legendre-Gauss-Lobatto 
points and weights in the x-direction and the z-direction. 

Solution of (!50|) . The solution x of (!50|) can be obtained by invert- 
ing the matrix {B^ — \fRC^). But this matrix has a large condition 
number. Thus we show a method to obtain the solution inverting the 
matrix given by f lT7|l which has a much smaller condition num- 
ber. For example, for A^^; = 10, Nz = 8, the condition number of 
(B^-VRC^) is 0(10^*^) while the condition number of is 0(10^). 

Since $i G -E'2 = span{0i, 02}"'"; we look for a solution of fl50|) in the 
form 

N 

(54) X = y^^XjXj, 

i=3 

where the eigenvectors of 

(55) B^Xi - y/RC^Xi = (5,[R)D^Xi. 

If we multiply dSO]) by [D^)"^ and use ([55]), the left hand side of (EOl) 
becomes 

N 

(56) ^a;,A(i?)a;, = (D^)"i(i?^ - v^C^)a; = [D^'r^ := /. 

We determine / from f = b using Gaussian elimination. Once 
again using Gaussian elimination, we can find the coefficients fi in the 
expansion 

TV 

(57) f = J2f^X^. 

1=1 

In f lFTj) . we see that /i = /2 = is necessary for the existence of 
a solution of ( 150|) . From (!56|) and ( 157|) . one finds Xj = fi/f3i{R), i = 
3,4, . . . , N. Thus the Jacobi expansion coefficients in (149|) of the center 
manifold are given by 



Numerical computation of in fl27j) . We approximate an by 
(58) af, = G,(0f,<l>f,0f). 

The integrands in Gs{4>i ,(^1) polynomials of degree at most 
(3A^2;+9, 3A''2+9). Thus to replace the integrals in (158|) . one needs again 
(f^x + 5, + 5) quadrature points and nodes in the numerical inner 
product. The other coefficients aij in ( l27ll are approximated similarly. 
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Remark 6.1. We observed that increasing and Nz above = 
10 + 2k and = 8 only changes af- in the seventh digit when the first 
critical mode which has k rolls and the second critical mode has k + 1 
rolls in their stream functions. 

7. Numerical Results and Discussion 

We computed coefficients of the reduced equations for various Pr 
values ranging from 0.1 to 10^ at the ffist three critical length scales 
and at the critical Rayleigh numbers which are given in Table HI 

As proved in Theorem I4.H the coefficients an and 024 are always 
negative. In our numerical calculations, we encountered that ai3 is 
also always negative. But the sign of 022 and the sign of D3 depends 
on L and Pr and are given in Table HI 

For the first critical length scale Lc = 1.5702, we found that 022 and 
D3 changes sign from positive to negative between 0.04 < Pr < 0.05 
and 0.14 < Pr < 0.15 respectively. Thus the transition is as described 
in Figure E]( a) for Pr < 0.04, as in Figure [S](b) for 0.05 < Pr < 0.14 
and as in Figure [S](c) for Pr > 0.15. Thus the mixed modes can be 
stable when Pr < 0.14 but only the pure modes are stable points of the 
attractor when Pr > 0.15. 

For the second critical length scale Lc = 2.6611, we always observed 
that a22 < 0. However, D3 changes sign between 0.05 < Pr < 0.06. 
Thus the transition is as described in Figure [5]^b) for Pr < 0.05 and as 
described in Figure [5]^c) for Pr > 0.06. In particular, the mixed modes 
can be stable when Pr < 0.05 but only the pure modes are stable steady 
states when Pr > 0.06. 

For higher critical length scales (third and beyond), we found that 
a22 < and < for the Prandtl numbers we considered. Thus the 
transition is as described in Figure [5](c). For this length scale, either 
the critical Prandtl number that was observed for the ffist two critical 
length scales is now very close to zero or it does not exist at all. 

The above analysis depends on the coefficients aij of the reduced 
equations and predicts the transitions when both eigenvalues P2 
are close to zero. Now we present an analysis depending on the di- 
rect computation of the numbers Di, D2 (both of which vanish when 
f^i = (^2 = 0) and D3. We computed Di, D2 and for L and 
R values around (but not necessarily very close to) the criticality 
(L, R) = {Lc, Rc) for the ffist three critical length scales and for Prandtl 
numbers Pr = 0.1, 0.71, 7, 130. The results are shown in Figure [6l Al- 
though we might have omitted the smallness assumptions of \L — Lc\ 
and \R — Rd where our main theorem is valid, these figures help us 
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Pr=0.1 



Pr =0.71 



Pr =1 



Pr=130 



3205.6 



lllll, 









f^65 1.570 1.575 1.565 1.570 1.575 1.565 1.570 1.575 1.565 1.570 1.575 
2207.8 




2:656 2.661 2.666 2.656 2.661 2.666 2.656 2.661 2.666 2.656 2.661 2.666 
1995.1 



19053!^ 



00 3.705 3.710 3.700 3.705 3.710 3.700 3.705 3.710 3.700 3.705 3.710 
♦ DUO, D2>0 • D1>0, D2<0 ■ D1<0, D2<0, D3<0| D1>0, D2>0. D3>0| 



Figure 6. The signs of Di, D2 and Da in the L-R 
plane. In each subfigure, the x and y axes denote the 
length scale L and the Rayleigh number R respectively. 
For each column, the Prandtl number is given above. 



predict the transitions in the L — R plane. The results we obtain are 
as follows. 

For Pr = 0.71, Pr = 7, Pr = 130, transitions are qualitatively same 
in the L-R plane. For L > Lc, the basic motionless state loses its 
stability to the eigenmode with wavenumber + 1 as the Rayleigh 
number crosses the first critical Rayleigh number and further increase 
of the Rayleigh number does not alter the stability of this steady state. 
This is in contrast to the situation L < Lc where there is a transi- 
tion of stabilities as the Rayleigh number is increased. Namely, as the 
Rayleigh number crosses the first critical Rayleigh number, the eigen- 
mode with wavenumber k becomes stable. As the Rayleigh number is 
further increased, both eigenmodes coexist as stable steady states and 
the initial conditions determine which one of these steady states will 
be realized. Finally as the Rayleigh number is further increased, the 
eigenmode with wavenumber k + 1 becomes stable. 

The transition at Pr = 0.1 is essentially different than for those 
at Pr = 0.71, 7, 130. In particular, for the first critical length scale 
Lc = 1.5702, for L < Lc, subsequently mode with wavenumber k, 
mixed modes and finally mode with wavenumber k + 1 will be realized 
as the Rayleigh number is increased while for L > Lc, k + 1 mode is 
the only stable steady state. 
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= 1.5702 


Lc = 2.6611 


= 3.7048 


Pr 


a-22 ^ iU 


Do V 1 




Do V 1 


Ploo V 1 (l2 


Do V 1 (1^ 


n ni 






9^? 79^19 


^ 4(1 98Q4 


9(1 ^il S^i 




n 04 


n 1 




fi (1784 


o.uuyu 


-o.ouyo 


1 441 


n 

yj .yjo 






-O.UO 1 ^ 


1 (14^1^1 




1 1 ^19^1 




-0 S740 

U.O 1 TiU 


90 09^,9 


-4 3Q09 


-0 971 3 


-3 8Q66 


-1 01 99 


14 


-2.3697 


0.6458 


-3.0486 


-2.3788 


-2.8749 


-0.9899 


0.15 


-2.4662 


-0.0707 


-3.0218 


-2.4414 


-2.8641 


-1.0182 


0.71 


-3.2469 


-9.5827 


-2.8940 


-2.5437 


-2.8848 


-1.4316 


7 


-0.8551 


-0.8348 


-0.7484 


-0.1610 


-0.7444 


-0.0987 


100 


-0.0687 


-0.0055 


-0.0599 


-0.0010 


-0.0595 


-0.0006 


1000 


-0.0069 


-0.0001 


-0.0060 


-0.00001 


-0.0060 


-0.000006 



Table 4. The coefficients of the reduced equations for 
various Pr values at the ffist three critical length scales. 
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